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Abstract 

We rigorously derive weakly nonlinear relation between cosmic density and velocity fields 
up to third order in perturbation theory The density field is described by the mass density 
contrast, 5. The velocity field is described by the variable 9 proportional to the velocity 
divergence, 9 = — f(Q)~ 1 H " 1 'V ■ v, where /(f2) — f2 ' 6 , Q is the cosmological density 
parameter and H is the Hubble constant. Our calculations show that mean 5 given 9 is 



a third order polynomial in 9, (S)\g = a\9 + ci2(# 2 — erf) + a^9 z . This result constitutes an 
extension of the formula (S)\g = 9 + a2(9 2 — erf ), found by Bernardeau (1992) which involved 
second order perturbative solutions. Third order perturbative corrections introduce the 
cubic term. They also, however, cause the coefficient a\ to depart from unity, in contrast 
with the linear theory prediction. We compute the values of the coefficients a p for scale-free 
power spectra, as well as for standard CDM, for Gaussian smoothing. The coefficients obey 
a hierarchy 03 <C ct2 *C a\, meaning that the perturbative series converges very fast. Their 
dependence on Q is expected to be very weak. The values of the coefficients for CDM 
spectrum are in qualitative agreement with the results of N-body simulations by Ganon et 
al. (1996). The results provide a method for breaking the f2-bias degeneracy in comparisons 
of cosmic density and velocity fields such as IRAS-POTENT. 

Key Words: cosmology: theory - galaxies: clustering - galaxies: formation - large-scale 
structure of the Universe 
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1 Introduction 



The most common assumption in theory of structure formation is the gravitational instabil- 
ity hypothesis: the observed large-scale structure has formed by the gravitational amplifica- 
tion of small-amplitude fluctuations present in the primordial density field. Cosmic velocity 
fields of galaxies result consequently from the gravitational attraction of large-scale mass 
inhomogeneities, that perturb the uniform Hubble flow. The quantitative relation between 
the peculiar velocity field, v, and the mass density contrast field, 5 = p/pt — 1, where pi, is 
the background density, can be inferred from the dynamical equations for the pressureless 
self-gravitating cosmic fluid. 

In linear regime, i.e. for 5 <C 1, the fluctuation field grows in time by an overall scale 
factor D(t) (which depends on the cosmological parameter Q), preserving its initial shape, 
<5(x, t) = D(t) 5(x,ti). As a result, the linear theory relation between the density and the 
velocity field is local 

5(x) = -/(fi)- 1 iJ - 1 V-v(x), (1) 

where H Q is the Hubble constant and f(Q) — fi ' 6 (see e.g. Peebles 1980). One can use the 
above formula to reconstruct from the large-scale velocity field the (linear) mass density 
field, up to an ^-dependent multiplicative factor f(Q). The comparison of the reconstructed 
mass field with the observed large-scale galaxy density field could therefore serve as a test 
for the gravitational instability hypothesis and as a method for estimating Q (Dekel et 
al. 1993). 

There are, however, both observational evidence and theoretical arguments for thinking 
that galaxies are biased tracers of mass. When the fluctuations are small one can assume 
that the galaxy and mass density contrast fields are linearly related, 5 g = b5, hence 

-^V-v(x) = ^(x). (2) 

The comparison of the POTENT-reconstructed mass field with the IRAS galaxy field yields 
/(f2)/&iRAs values close to unity (Dekel et al. 1993, Dekel 1994). It is then tempting to 
conclude that the large-scale dynamics is consistent with an assumption of Q = 1, provided 
that 6pRAs is a ls° close to unity. However, since we do not know anything a priori about 
bias, we should measure it independently. 

It has been suggested that nonlinear corrections to the linear density-velocity relation 
(hereafter DVR), equation (JJ), can help to perform such a measurement (Yahil 1991). The 
corrections are indeed necessary because there are points in the ^potent-^iras correlation 
diagram for which the density contrast reaches unity, clearly contradicting the underlying 
assumption of 8 <C 1. On the other hand the rms fluctuation of the mass field, a, is smaller 
(but not much smaller) than unity that means that the field is weakly nonlinear. 

In weakly nonlinear regime perturbation theory can be efficiently applied, as the re- 
sults of N-body simulations show (Juszkiewicz et al. 1995, Bernardeau 1994a,b, Baugh, 



2 



Gaztanaga & Efstathiou 1995, Lokas et al. 1995, Bernardeau & van de Weygaert 1996). 
However, most of the attempts to derive a weakly nonlinear extension of the linear DVR 
have been based on the Zel'dovich approximation and its modifications (Nusser et al. 1991, 
Gramman 1993). The Zel'dovich approximation is a useful qualitative guess of nonlinear 
dynamics but it provides only approximate answers for rigorously derived higher-order per- 
turbative solutions. Consequently, it does much better than linear theory, but still does not 
predict accurately the weakly nonlinear relation between velocity and density, as verified 
by N-body simulations (Mancinelli et al. 1994, Ganon et al. 1996). 

The first attempt, and so far the only one, to calculate DVR within the framework of 
rigorous Eulerian perturbation theory has been taken up by Bernardeau (1992a; hereafter 
B92). B92 has calculated the exact DVR for an unsmoothed final field in the limiting case 
of vanishing variance. The assumption a — > greatly simplifies mathematical calculations. 
It is not, however, well suited for the application to the IRAS-POTENT comparison: we 
are then not interested in the statistics of very rare events (5 ^> er) of linear field (a <C 1), 
but in the statistics of 'typical' events (5 ~ a) of a weakly nonlinear field (cr < 1). N-body 
cosmological simulations show that the exact formula of B92, when straightly applied to 
the case a < 1, works worse than the Zel'dovich approximation (Mancinelli et al. 1994, 
Ganon et al. 1996). 

B92 has also computed the first nonlinear (i.e. quadratic) correction for the DVR in the 
case of a smoothed final field with non-vanishing variance. However, neither the details 
of the derivation, nor the explicit form of the coefficient of the corrective term are given 
in the paper. On the other hand, a perturbation theory-inspired approximation of density 
as a third-order polynomial of velocity divergence turns out to be an excellent robust fit 
to N-body results (Mancinelli et al. 1994, Ganon et al. 1996). Theoretical construction of 
such a polynomial requires third order perturbative solutions and provides therefore higher 
order corrections to the DVR than those given by B92. All this inspired us to calculate the 
weakly nonlinear DVR of a smoothed final field with a ^ 1 up to third order in perturbation 
theory. 

The paper is organized as follows: in section 2 we derive weakly nonlinear DVR in its 
general form. In section 3 we compute values of the coefficients entering this relation for the 
case of scale-free power spectra, as well as for standard CDM. Discussion and concluding 
remarks are given in section 4. 

2 General derivation of the density-velocity relation 

In perturbation theory, one expands the solution for the density contrast as a series around 
the background value 5 = 0, 

S = <$! + 8 2 + S 3 + . . . , (3) 
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and truncates it at some order. The linear theory solution mentioned in section 1 is just 
the perturbation theory series truncated at the lowest, i.e. first order term, 



6 1 (x,t) = D(t)6(x,t i ). (4) 

Higher-order solutions are found iteratively: the second order contribution 62 is the solu- 
tion to the dynamical equations with 5\ as the source term for nonlinearities, and so on. 
Throughout this paper, we will consider only the growing modes, as the remaining ones 
are suppressed during linear evolution. In general, the n-th order solution is found to be of 
the order of (5i) n (Fry 1984; Goroff et al. 1986). Let us define a as the square root of the 
variance of the linear density field, i.e. a 2 = (5f), with (•) meaning the ensemble averaging. 
We have 5\ ~ a, 5 n ~ a n , and the series @ is a power series in a small parameter a. 
We describe the velocity field by a variable proportional to the velocity divergence, 

9(x,t) = -f(n)- 1 Hu 1 v-v(x,t) (5) 

(which is slightly different from the commonly used definition, e.g. Bernardeau 1994a). The 
variable 9 is as well expanded in a series 

= 1 + 9 2 + e 3 + ... . (6) 
The linear theory solution, equation (JJ), therefore gives 

*i(x) = *i(x). (7) 

Second order contributions to 5 and 9 are different. Their explicit dependence on Q 
is extremely weak and in the range of cosmological interest, 0.1 < Q < 3, the solutions 
are excellently approximated by the expressions which hold in the case of the Einstein-de 
Sitter universe (Bouchet et al. 1992), namely (Goroff et al. 1986) 

5 2 (x, t) = ^ 5? + dj, d a A x + ^d a dp A^fyAi , (8) 

and 

2 (x, t) = j 51 + dj x a Q A x + jd a d p AAfyAx . (9) 
Here, Ai(x, t) is the linear gravitational potential, 

We see therefore that up to second order in perturbation theory the divergence of the 
velocity field, V • v(x) = —f(Q)H (9i + 9 2 + ■ ■ ■), depends explicitly on Q only via a 
multiplicative factor f(Q). 

Equation (^) means that if we read from a linear field the values of pairs (5(x), 9(x)), 
point by point, and plot them on the 5-9 plane, they will lie on a straight line. The second 
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order contributions to <5(x) and 9(x), in addition to the local term ~ 5i 2 (x), contain non- 
local terms due to the linear gravitational potential, ~ QVx"^i(x')5i(x"). As a result, 
the weakly nonlinear DVR is no longer local. However, given 9, the spread in the values of 
S comes only from nonlinear corrections. Consequently, the points on the 5-9 plane are still 
strongly correlated: they are expected to form an elongated set of length ~ a and width 
~ a 2 around the mean trend (B92). This has been also observed in N-body simulations 
(Nusser et al. 1991, Bernardeau & van de Weygaert 1996). The mean trend can therefore 
serve as a very useful local approximation of a true nonlocal DVR. 

Full information about the density- velocity correlation is contained in a joint probability 
distribution function (PDF) of weakly nonlinear variables 5 and 9. Pioneering works on 
computing PDFs of weakly nonlinear cosmological density and velocity fields have been 
performed by Bernardeau. First, Bernardeau (1992b) computed a one-point PDF of an 
unsmoothed weakly nonlinear density field. Subsequently, he extended his calculations 
for the case of a top-hat window function (Bernardeau 1994a), and computed as well the 
PDF of a top-hat smoothed velocity divergence field (Bernardeau 1994b). A joint density- 
velocity PDF, however, is still known only for the case of an unsmoothed final field, and in 
the limit o — > 0. B92 calculated it in the form of mean 9 given 5; this relation is however 
easily invertible and the result is 

/ 2 \ 3/2 

5=[l + -9) -1. (11) 

Note that in the linear theory limit, 9 <C 1, the above equation indeed reduces to equa- 
tion ©. 

Juszkiewicz et al. (1995) by means of N-body simulations have shown that a one-point 
PDF of a single variable 5 (or 9) in the range of 'typical' events 5 ~ a can be very well 
approximated by the so-called Edgeworth series (e.g. Longuet-Higgins 1963, 1964 and refer- 
ences therein). The Edgeworth series is constructed from cumulants of the true distribution, 
defined as the connected part of the moments, 

^n \P ) conn • (12) 

The cumulants of order n > 2 provide an effective measure of non-Gaussianity because 
for a Gaussian distribution they vanish. Throughout this paper we will assume Gaussian 
initial conditions. Consequently, all n > 2 cumulants of a fluctuation field are initially 
zero. During nonlinear phase of evolution, however, they acquire nonzero values. Fry 
(1984) showed that cumulants of cosmic density and velocity fields in weakly nonlinear 
regime obey the following scaling 

K n = S n a 2 ^ + 0(a 2n ) . (13) 

To calculate the coefficient S n , the perturbative solution of (n — l)th order is needed. Let 
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us define dimensionless quantities related to cumulants as follows 



K 



= ~^ , (14) 



n/2 
2 



where by definition k 2 = (S 2 ) = cr| is the nonlinear variance of a density field. From 
equation (|T3|) one can deduce the order of weakly nonlinear corrections to its linear value 

a 2 = a 2 + G(a 4 ) . (15) 

The coefficients A„ are thus the cumulants of a standardized variable \x = S/us and we will 
refer to them as to 'standard cumulants'. The first two nontrivial standard cumulants, A3 
and A4, are called in statistics skewness and kurtosis, respectively. 
The Edgeworth series reads 



V Z7T 



1 + -\ 3 H 3 (fi) + ^A 4 # 4 (/i) +^A 2 ^ 6 (/i) + 
where if n (/i)'s are the n-th order Hermite polynomials generated by 



(16) 



(-l) n ^e-^ 2 = e-fHM . (17) 

In Table 1 we provide explicit forms of the few lowest order polynomials. From Equa- 
tions flTJ, (UJ) and (|TJ) we have 

A„ = S n a n s - 2 + 0(a 5 n ) , (18) 

which expresses the scaling behaviour of standard cumulants of a weakly nonlinear density 
field evolving from Gaussian initial conditions. In particular, A3 = S^as and A4 = S^a 2 , 
i.e. during weakly nonlinear evolution skewness and kurtosis grow like the rms fluctuation 
of the field and the square of it, respectively. In cosmology, there is a long tradition to 
call 'skewness' and 'kurtosis', respectively, the coefficients 5*3 and 5*4 themselves. We will 
honour it hereafter. 

Equation ( |I8|) ensures that the Edgeworth series is a series expansion of an exact PDF in 
powers of a small parameter as (Longuet-Higgins 1963). In weakly nonlinear regime we can 
thus approximate the true PDF by the Edgeworth series truncated at some order. Using 
equation ( 18|) the Edgeworth expansion, equation fllBD , can be rewritten in the explicitly 



perturbative, third-order form 

1 1 1 _ . 1 

(19) 



p(aO = -75= ^ 12 



1 + l -Sza s HM + Ls^lH^) +lslc7 2 s HM 

The Edgeworth expansion for the variable 6 or v = 9/ae, where a 2 = (8 2 ) is the nonlinear 
variance of the velocity divergence field has the same form, except that S3 and S4 are then 
the skewness and the kurtosis of the velocity divergence field. 
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The third-order Edgeworth expansion describes accurately the shape of a true PDF up 
to /i ~ aj l (Juszkiewicz et al. 1995). The failure of the approximation in the very tails 
reflects the fact that it is constructed only from a few low order cumulants of the true 
distribution (see also Bernardeau & Kofman 1995). In the present paper, however, we are 
not interested in the statistics of very rare events in the 5-9 space. Instead, we want to 
calculate just the lowest conditional moment: mean 5 given 9. For this purpose, we have 
to know the approximate form of the joint distribution for the variables 5 and 9 that needs 
to be accurate only for typical events, 5 ~ as, 9 ~ a$. The two-point generalization of the 
above third-order Edgeworth series exactly satisfies this condition. 

In fact, one can proceed in two ways. One can derive joint Edgeworth expansion and 
then calculate conditional moments from it. One can also, however, calculate the moments 
directly. Deriving the third-order joint Edgeworth expansion is a straightforward, but 
lengthy calculation, while of most interest for cosmology is just the first moment, describing 
the mean trend. Therefore, in this paper we calculate it directly, postponing the calculation 
of the joint Edgeworth expansion and higher-order moments (e.g. the variance around the 
mean trend) resulting from it to the next paper. The methods of calculating moments and 
the full PDF are still closely related and in the following calculation we are inspired to some 
extent by Longuet-Higgins (1963), who derived a second-order joint Edgeworth expansion 
in order to apply it to statistical theory of sea waves. 

The conditional probability for 5 given 9 is 

(20) 

where p(d~, 9) is the joint PDF for 5 and 9. The characteristic function of p(5, 9) is 

$(zt,zs) = J J e itS+lse p(S, 9)d5d9 . (21) 
Expanding the exponentials we obtain 



*(it,is)= £ W(^) m («)"> (22) 

m,n=0 m]n] 

where (5 m 9 n ) are the joint moments of 5 and 9, 

(S m 9 n ) = J J 5 m 9 n p{5, 9) d5&9 . (23) 
If the joint moments are known, p(5, 9) can be calculated via the inverse Fourier transform, 

p(5, 9) = —- J I e- it5 ~ tsd $(it, is) dtds . (24) 
Mean 5 given 9, (o)\d, is by definition / 5p(5)\gd5. From equation (|20|) we have 



S5p(5,9)6S _ 
{6)le= P (9) ■ (25) 
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Let us denote / Sp(8, 9)dS by Af. By equation ( p4|) 

1 



(27T) 



-it5—isd 



5® (it, is) dtdsdd 




(2tt 
1 

2^ 



7T) 2 7 7 a(«t) j 



d 



e- ts "$(it, is) -^-S D (t) dtds . 



(26) 



where 5^(t) denotes the Dirac delta function. Integrating by parts we obtain 



Af = — 
2vr 



-isO 



d 



d(it) 



§(it,is)\ t= o ds. 



(27) 



The characteristic function is related to the cumulant generating function, JC, by the 
equation 

$(it,is) =exp[}C(it,is)]. (28) 



The cumulants, n mn , from which /C is constructed, 



K 



E 



K , 



mini 



:(it) m (is) 



(m,n)^(0,0) 

are given by the connected part of the joint moments 



(29) 



(5 m e n ) c 



(30) 



Using equations (p8|) and (|29|) we obtain 



<9 



$(it,is)\ t=0 



.71=0 ' L - 



exp 



' 00 

.n=l 



(31) 



By definition, «on are the ordinary cumulants of the variable 9. The variables S and 9 have 
zero mean, so k w = k 01 = 0. Equations (|27| ) and (|3~TD then give 



AT 



2vr J- 



ds e 



-isB 



.n=l 



exp 



E 

,n=2 



US 



1 rOG 

2tt J -00 



.n=l 



exp 



00 

E 

.n=3 ' 4 - 



(32) 
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Let us define a new variable z = k 02 s and let us recall that \i and v are the standardized 
variables, 
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We then have 



Af 



2tik 



1/2 
02 



°8 ^ 



and v 



9 9 



/t 1/2 ' 

^ ^02 



E -^j, («r 

.71= 1 TIIKQ2 



exp 



E 



^ , n/2 

.n=3 n!« 2 



(33) 



(34) 



The standard joint cumulants are defined by 



ft i. 



m/2 n/2 ' 



(35) 



K 



20 K 02 



hence, K 0n /n% 2 



\ l n / 2 

AQn, K ln/ K 02 



U = ± (!*) 



1/2 



= «2o 2 -^in and 

dze ~^+2iuz) 



00 \ 



.n=l 



ri! 



exp 



00 \ 
—(iz) 



n=3 



n ! 



(36) 



One may ask why we have introduced the cumulant generating function: using just the 
characteristic function $, the above equation would look formally simpler. The reason is 
similar to that in case of constructing one-point PDF. From perturbation theory it follows 
that standard joint cumulants, equation fl35"D, obey the following scaling hierarchy 



A, 



C m+n— 2 1 fs} ( ^m+n \ 



+ 0( 



a 



(37) 



where a is the linear variance of 5 or, equivalently, of 9 (recall that at linear order 5 = 9). 
The series in equation (|36D are therefore power series in a small parameter a and truncating 
them at some order p we neglect contributions which are ~ a p+1 . Perturbation theory also 
predicts that 

°l = (0 2 ) = k 02 = a 2 + 0{a*) , (38) 

so when we are interested in the leading order terms in hierarchy fl37|) we can use linear o 
instead of nonlinear a e (or ag, see eq. ]I5| ). 

In the present paper we want to calculate the weakly nonlinear extension of the linear 
5-9 relation, up to cubic in 9, 0(<Tq) terms. In equation (|36D , relating mean /1 = 5/as and 
v = 9/<jg, we will thus keep terms up to the order of of. We have 



dz e 



An^) + — iiz) + — 



-Uz 2 +2ivz) 



IZ) 



1 , ^03/. \3 



A 



^1(^)4 + ^03(^)6 

24 v 1 72 K J 



. (39) 



In the expression above, A 



.12 



A 



03 



ao and A 



13 



A 



01 



\ 2 '2 

A 3 ~ a e 



The cumulant An 



deserves a separate, more detailed treatment. Defined as (59) /[{5 2 ) l/2 {0 2 ) 1 ' 2 ] [see eq.|po| J, 
it is the correlation coefficient between the fields 5 and 9. Since the fields are identical 
at first order, 5\ = 9±, at the lowest order An = 1. From equation ( p7|) it follows that 
the higher-order correction to this value of An is 0(ctq), so in general An = 1 + 0(aj). 
Multiplying the polynomials in equation ( p9[ ) we keep only the terms up to the order of 
of. It means also that we replace the products AnA mn with m + n > 3 by A mn , since the 
correction is of at least cubic order in <jg. After sorting the resulting terms of the form 
(iz) n we integrate them, using the identity 

1 



'27T 



dze~i {z +2ivz) (iz) n 



HJv) 



1,-2 
2' 



(40) 



9 



where H n are the n-th order Hermite polynomials. The result is 

1 /^o\ 1/2 



e> 2 



x 



+ 



fhx \K 02 J 

2 o 



A 



13 



6 



A()4 

24 



AitA 



12^03 
12 



X 2 
72 



(41) 



To calculate mean S given # from equation ( [25] ) we need also one-point PDF of velocity 
divergence. As already stated, it is given by the one-point Edgeworth series, equation fllCf), 



1/2 



for the variable v = 9/ k 02 
p(v) -- 



^7T 



1 + ^A 3# 3 (^) + ^A 04 ^4(^) +^\ 2 03 H 6 (v) 



(42) 



Recalling that = S/k^q 2 , and by equation (p5|), we have (/x) | ^ = ^2o 1 ^ 2 ( ( ^)le = K 2o ^ N /p(9) 
= (n 02 / k^Y^M /p^) . From equations (f4l|) and ([42]), expanding the denominator, multi- 
plying, and keeping only the terms up to C(o"|) we finally obtain 



where 



^03 



V 2 {v) = > fH 2 {p)+ f) 



(43) 
(44) 



and 



6 3V ; 24 



+ 



X 2 
A 03 

72 



# 7 M 



z/if 6 (i/) - 2H 3 (u)HJv 



AitA 



12^03 
12 



H 5 {v) - H 2 {v)H 3 (v) 



2uH 2 3 (u) 



(45) 



Note that T^^) is and V 3 (v) is 0{oq). (Remember that v itself is standardized, so 

v ~ = C(l) and erg-dependence comes only from the standard cumulants.) 

Equation fl4"3] ) expresses weakly nonlinear density- velocity relation (DVR). In linear re- 
gime, oq — > 0, the correlation coefficient between 5 and 9 is unity, An = 1. Moreover, 
in this limit V 2 and V3 also approach zero, so from equation (]43| ) we reobtain the linear 
theory result, (fj)\ v = v, or = 9. The polynomials P2 and V3 are higher-order 

corrections to this linear relation. As we took into account the corrections up to third 
order in perturbation theory, we do not expect terms of order higher than cubic in v to 
appear in equations (|44])-(|45]). Indeed, in these equations all the terms of the form v n with 
n > 3 remarkably cancel out. To prove this, we use the recurrent relation for the Hermite 
polynomials 

H n {u) - uH^iy) = -(n - l)H n „ 2 {v) , (46) 
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along with their explicit forms for n — 1, . . . , 5 (see Table 1). The result is 



Al2 — Aq3 / 2 



and 



VM = Al3 : Ao4 (z/ 3 - 3u) + (Ao3 ~ Al2)A ° 3 (z/ 3 - 2;/) 



(47) 



(4£ 



6 2 
We see that V^v) and ^(f) contain quadratic and cubic corrections in u, respectively. 

Let us now deal with the joint cumulants in V%. We have \qzo\ = k 03 = (9 3 ) = Ssqctq, 
where S30 denotes the skewness of the variable 9. Thus we get 



Aq3 — S30OQ . 



(49) 



The other, mixed cumulant in Viiy) is defined as Ai 2 <7g = ((S\ + 5 2 + ■ ■ -)(9i + 9 2 + ■ ■ -) 2 ). 
Recalling that 8\ =61, we obtain 



A 



12 



<5Q 



(50) 



where S%s denotes the skewness of the variable 5. Thus A12 is a linear combination of 
ordinary third-order cumulants of the single variables 5 and 8. Equations fl4"T|) , fl49| ) and (p0|) 
then yield 



AS, 



6 



<x^ 2 - 1) , 



where 



AS* = S 



38 — S38 . 



(51) 



(52) 



We can now calculate the lowest order weakly nonlinear extension of the linear DVR. 
As we will show later, An = 1 + 0(oq). Keeping the terms up to 0(oq) in equation ([43|) 
we thus have 



(/i)|„ = l + 0(al 



AS 2 
~6~ 



a e (v 2 -l)+0(al 



or 



AS, 



(li)\ v =v + ^a 6 (p*-l)+0(al). 
o 



(53) 



(54) 



^From equation (|38|) it follows that os = o~e + 0(ag), so fi = 5/<js = + 0(&g), hence 

A 9 

(6)1 = a e u + —^°l(v 2 - 1) + O(ol) . (55) 
o 

Since v = 9/a e we end up with 



(S)\ e = 8 + 



6 



0(al 



(56) 



The above equation is the lowest, second order weakly nonlinear extension of the linear 
DVR. Consequently, the coefficient of the corrective term is composed from cumulants 
calculable at second order (535 and S 3 g), and the term is quadratic in 9. This term is shifted 
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down additionally by aj. Note from equation (|25|) that J (5)\gp(9)d8 = J 5 p(5,9) d5d9 = 
J 5p(5) d5 = (5) (not conditional, but ordinary), that is zero by definition. The <r| term, 
naturally emerging from our calculations, precisely ensures this. 

Apart from deriving the exact DVR in the case of an unsmoothed field with vanishing 



variance, equation (|TT|) , B92 calculated also the second-order DVR including the effects of 
finite variance and smoothing of a final field. Our result, equation (pB[), coincides exactly 
with equation (17) of B92 with the coefficient B = AS 3 /6. 

We will deal now with the cumulants in the V 3 term. From equation (f49|)-(|50"l) we have 

(A03 — Ai2)Ap3 _ AS 3 S30 2 /c-yN 

2 ~ 6 ° B ' [ ' 

The joint cumulant A13, unlike A12, is not a linear combination of ordinary cumulants of 
the single variables 5 and 9. We have (A 13 - A 04 )^ = ($9 3 ) - (9 4 ) = ((5 - 9)9 3 ) = 
({6 2 - 9 2 + 5 3 - 9 3 + • • -)(0i + 9 2 + 9 3 + • • -) 3 > and therefore 

A13 — A04 = S 4 o"g , (58) 

where 

E4 = m02(5 2 - 9 2 )) + (5%) - (9f9 3 ) ^ 

cr 6 

In the expression above (■) stands for the connected part of the moments. Note that E4 is 
not equal to S45 — S40: while the last two terms in equation ( p9|) are indeed parts of the 
expressions for the ordinary kurtosis of a single variable 5 or 9, respectively, the first term 
is a truly mixed moment and constitutes a new quantity. Using the results d57|)-(|58D in 
equation (pE8|) we obtain 



£ 4 — AS 3 S 39 2 3 



AS3 S 3 e _ S4 
3 2 



ajv . (60) 



Equation fl43|) expresses weakly nonlinear extension of the linear DVR up to Oiaf) 
corrections. The scaling of the standard cumulants with o~ 9 , equation (|37D, ensures that 
it was enough to calculate A mn with m + n > 3 at the lowest order. The corrections to 
An, however, are O(o"|), so they cannot be neglected. Similarly, fi = 5/k 1 Jq = 5 /as = 
(cg/as)(S/a0) = [1 + 0(al)]5/ag, so the corrections to the linear evolution of the variance 
of 6 and 9 should be taken into account as well. Changing the variables in equation (fl3|) 
to 5 and 9, equation (|3~3D, we have 

k- \ 1/2 



(S)\e=( — ) Xn9 + ci0[V 2 (9/a0)+V 3 (9/a0)} + O(al). (61) 

^02- 



By definitions (0) and (BO), 



^20 \ ^ . «u _ (69) _ (59) - (9*) + (9*) _ ((5-9)9) 
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which after expanding 5 and 9 in perturbative series gives 



^02 



A u 0=[l + £ 2 ^]0, (63) 



where 

^ = - ft)) + (W - (W ^ (g4) 



a 1 

f2\ //)2\ 



Again, E 2 is not equal to (5 ) — (9 ). The last two terms in equation ( pi|) are indeed parts 
of the expressions for nonlinear corrections to the linear evolution of variance of S and 9, 
but the first term is a truly mixed moment and constitutes a new quantity. 



To obtain third-order weakly nonlinear DVR in its final form we combine equations (|5lD 
(^) and ( |B"3"D with fl5"T|) ■ Note that "P 3 contains also a term linear in 9. We have 



(5) \ e = a x Q + a 2 (# 2 - a 2 e ) + a 3 9 3 , (65) 

where 



ai = 1 + 



AS3 5 3 e £ 4 

ZjO H — 

3 2 



(66) 



a 2 = (67) 
6 

as = — ^~ LJi , 68 

6 

with £ 2 and E 4 given by equations (0) and (|59|), respectively Equations (165])- (|68|) con- 
stitute the main result of this section. Note that we reobtain the second-order DVR, 
equation ([56]), when we neglect the 0(cr|) terms (i.e. ~ a 2 9 and ~ 9 3 ). 



An important conclusion can be drawn immediately: the DVR of a weakly nonlinear 
{o~e 1) field is different from the linear theory prediction, equation (0), even for \9\ <C 1. 
Namely, in this case 

(S)\ e = ai 9-a 2 a 2 e . (69) 

Thus, the linear relation is shifted down, as already discussed. What is perhaps even more 
interesting, the coefficient a\ generally departs from unity. The strength of this shift and 
departure depends however on the particular values of the coefficients a n . This is the 
subject of the next section. 



3 Numerical calculations 

A brief outline of the perturbative solutions to the Newtonian equations of motion needed 
for numerical calculations of the coefficients a n is given in Appendix A. 

The smoothing of the fields on scale R is introduced by the convolution of the density 
contrast (or velocity divergence) field and the filtering function W 

6 R (x,t) = Jd 3 y8(y,t)W(\x-y\,R). (70) 
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We perform our calculations for a Gaussian filter function which is appropriate for obser- 
vational analysis of cosmic velocity fields and comparing them with the density fields (e.g. 
POTENT and IRAS-POTENT comparison). The Fourier representation of the Gaussian 
window function is given by 

W(kR) = e~ k2R2/2 . (71) 
We assume a Gaussian distribution for the first order 5i and 61 and define 

<x 2 = (51) = D\t) J P(k) W\kR) (72) 

as the linear variance of the density (velocity divergence) field. We assume that for a < 1, 
the first few terms in the perturbative expansion provide a good approximation of the exact 
solution. Since S± and 9\ are assumed to be Gaussian random fields, all their statistical 
properties as well as those of the higher order terms in the perturbative series are determined 
by the power spectrum P(k), defined as 

(Mp)Mq)> = (2n) 3 5 D (p + q)P(p). (73) 



3.1 The calculation of the coefficient 

The values of skewness for density contrast and velocity divergence fields, given to the 
lowest perturbative order respectively by 

a, = 5#> (74) 

(j 

and 

S„ = ^, (75) 

depend on the assumed form of the power spectrum. 

We begin by considering spectra with a power-law form 

P{k) = Ck n , -3 < n < 1, (76) 

where C is a normalization constant. For such fields, smoothed with a Gaussian filter, the 
linear order contribution to the variance given by equation QT2j ) is 

° 2 = cd2 VjS^' (77) 

where R is the smoothing scale. The values of skewness are (Lokas et al. 1995) 

„ fn + 3 n + 3 3 1\ / 8\ „ fn + 3 n + 3 5 1\ 
Su = 3 ,F, (— — -, j) - („ + -) 2 F, (— , — , -, j) (78) 

fn + 3 n + 3 3 1\ / 16\ „ fn + 3 n + 3 5 1\ 
S M = 3 2 F 1 (— — )_(„+) ^ (—,— -,-) (79) 
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where 2 -Pi is the hypergeometric function. Therefore the coefficient a 2 is 

S35 -Sm 4 fn + 3 n + 3 5 1\ 

The result is very weakly dependent on the value of the Q parameter (for details see the 
appendix in Lokas et al. 1995; Bernardeau et al. 1995 and Bouchet et al. 1992). A good 
approximation for the Q dependence in the range 0.1 < Q < 3 is obtained by replacing the 



constant coefficient 4/21 in equation fl80 ) with an ^-dependent function 



G(n) = i-xspjm (S1) 

where 

K(Q) = ^- 2/63 , C{Q) = ^Q- 1/21 . (82) 

The second column of Table 2 gives the values of 02 for integer and half-integer values of 
the spectral index n and = 1 while Figure 1 shows the coefficient 02 as a function of n 
for three different values of fl 

We have chosen the scale-free spectra of the form ([76|) not only because of their simplicity 
but also for their straightforward applicability to realistic power spectra. Indeed, in the 
case of higher order cumulants the value of the cumulant (the skewness or the kurtosis) 
is very well approximated by the result for the scale-free spectra with the effective index 
defined as (Bernardeau 1994a) 



_ R do*(R) 

neff -~^^R-- 3 - (83) 



As an example of a scale-dependent power spectrum we consider the standard CDM 
spectrum 

P(k) = — w (84) 

{1 + [hk/r + {hk/vf 2 + {hk/v) 2 ] v } lv 

with n — 1, r = 0.5, v = 1.13 and the constants in units of h~ x Mpc are l\ = 6.4, I2 = 3.0, 
I3 = 1.7 (e.g. Efstathiou, Bond & White, 1992). We normalize the spectrum so that the 
linear rms density fluctuation in spheres of radius R = 8h _1 Mpc is equal to unity. Thus 
the definition of variance ( ff^) together with the following shape of the spherical top hat 
window function in Fourier space 

W TH {kR) = ^{kR)^' 2 J V2 {kR) (85) 



(where J is the Bessel function) and the power spectrum ( [84]) yield the normalization 
constant of C — 4.09 x 10 5 {h~ x Mpc) 4 . Note, that the normalization procedure is the only 
place where we use the top hat filter, all other calculations are performed for a Gaussian 
window function (ffll). 
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We have calculated the coefficient a 2 for the CDM spectrum at different smoothing 
scales R in two ways. First we have found numerically the exact values of Ss$ and S30 for 
CDM spectrum following the procedure of Lokas et al. (1995). These values are shown in 
the third and fourth column of Table 3. They produce the exact value of the coefficient 
a 2 given in the fifth column of the Table. The second option was to use the formula (jS0|) 
with n e ff corresponding to each scale, calculated from equation fl8"3"|). The second column 
of Table 3 gives the values of the effective index corresponding to each of the scales. Thus 
obtained value of a 2 is given in the last column of Table 3. The comparison between the 
two values of ci2 shows clearly that the discrepancy between them is less than 1% at all 
scales. The exact values of a 2 for the CDM are repeated in Table 4, which summarizes the 
results for this spectrum. 



3.2 The calculation of the coefficient 

As we have shown in the previous section, the coefficient 03 is given by 

as = — 86 

6 

where 

S4 = H$s 2 e 2 ) - me*) + - <m) _ (87) 

We have named the quantity £4 to stress its similarity to the kurtosis of density and velocity 
divergence fields, which to lowest order in perturbation theory are given respectively by 

and 

6<gg) + 4(fffl 3 ) fM . 
£>4e = (»yj 

(unless the initial conditions are non-Gaussian: for details in the latter case see Chodorowski 
& Bouchet 1996). This shows that most of the expressions constituting the value of S 4 for 
power law spectra has already been calculated by Lokas et al. (1995) while performing the 
calculations for kurtosis. Since they were not published, we give them in Table 5 for integer 
and half integer values of the spectral index n. They will also be needed in the calculations 
of the next subsection. The only unknown part of £4 is the expression of the form ( 618262) 
which is calculated in Appendix B. 

In the case of no smoothing (when window function W(kR) = 1), which also corresponds 
to putting the spectral index n = —3, a completely analytic result can be obtained fairly 
easily; we get (8f8 2 6 2 )/a 6 = 1768/441 ps 4.01. In this case we have 

"• = -s^ a0101 (90) 
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The sixth column of Table 5 gives the values of S 4 calculated according to (|87|) . Finally, 
the third column of Table 2 lists the values of 03 obtained from equation (|8~6D in the case 
of power-law spectra for integer and half integer values of the spectral index n. 

Bernardeau (1994a) showed that third order solutions for 5 and 8, similarly to second 
order, depend very weakly on Q. Since the coefficient (13 is constructed from terms up 
to third order, its expected dependence on Q is very weak. Bernardeau (1994a) did not 
give explicit forms for weakly ^-dependent third order solutions. Unlike to the case of 
the coefficient a 2 (see previous subsection), we cannot therefore verify in detail the above 
supposition. Still, we are able to do this at least for the case of the spectral index n = —3. 
The exact formula of B92, equation (|TT|) , describes the limiting case a 2 — > 0, n = —3 and 
£1 = 0. One cannot thus use it to deduce the value of the coefficient a\ which includes 
corrections 0(<j 2 ). On the other hand, one can use this formula to derive the values of the 
coefficients a 2 and a 3 , which are calculated in the limit a 2 —>■ 0. The Taylor expansion of 
equation (|TT|) yields 

a 2 (Q = 0) = -« 0.167 (91) 
6 

and 

a 3 (Q = 0) = -— w -0.0185 . (92) 
54 

The corresponding values for a 2 and a 3 in the case n = —3, Q — 1 are respectively 0.190 
and —0.0101 (Table 2). The relative change of the value of the coefficient 03 is therefore 
greater than the relative change of a 2 . Nevertheless, also 03 depends on Q extremely weakly 
in a sense that it almost vanishes both for £7 = 1 and £7 = 0. All kurtosis-type quantities 
entering the definition of 0,3 are of order of unity (see Table 5) and very precise cancellation 
of them is needed to assure (13 -C 1. Therefore even weak dependence of the perturbative 
solutions on Q could in principle destroy this 'fine tuning'. For n = —3 this is clearly not 
the case. 

In fact we were able to check it rigorously for all values of n for the kurtosis-type 
quantities in E4 (eq. |87|]) that involve only second order solutions. In the range of 0.1 < Q < 
3 we found that the ^-dependence of (S 2 5 2 9 2 ) and {6\8 2 ) is similar and almost cancels out 
when the two are subtracted. Analogously, Bernardeau (1994a) noted that the combination 
S40/S30 is almost independent on Q, to much bigger extent than the moments £30 and S40 
themselves. Very weak dependence of 03 on Q is an interesting problem and we will address 
it in more detail elsewhere. 

To obtain the values of the coefficient a 3 for the CDM spectrum we apply the effective 
index method described in the previous subsection. For each of the indices calculated from 
equation (p3[) for a given Gaussian smoothing scale we interpolate the value of a 3 from the 
values given in Table 2 using an accurate polynomial fit. The results are presented in the 
last column of Table 4. 
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3.3 The calculation of the coefficient a\ 

It has been proved that the coefficient a\ is given by 



ai = 1 + 
where 



2 



< (93) 



Ss = <w - m) + (w - (w (94) 

and £4 has been defined and discussed in the previous subsection. 

The quantities involved in £ 2 are of the form of the lowest order weakly nonlinear 
corrections to the variance of the density and velocity divergence fields which are of the 
order of a 4 (Lokas et al. 1996) 

°l~° 2 _ <^ 2 > + 2<W (95) 



a 4 a 4 



a 2 - a 2 (0f) + 2{<khl (96) 



a 4 a 4 



We recall that a 2 and <jg stand for nonlinear variance of density and velocity divergence 
respectively while a 2 is the linear variance given by equation ([721). In equation ( |9~3"D a 2 can 
be replaced by a 2 since their difference is already O{of). 

The details of calculations of the terms involved in £2 are given in Appendix C. As 
discussed in the Appendix, some of the terms are divergent at spectral indices n > — 1. 
Instead of dwelling on those divergences we focused our attention on analysis concerning 
scale- free power spectra in the case of no smoothing and the case of —2 < n < — 1 with 
smoothing, which is well justified observationally. As recent analyses of measurements 
suggest (Gaztanaga 1994; Feldman, Kaiser & Peacock 1994; Peacock & Dodds 1994) the 
linear power spectrum can be approximated over large range of scales by a power law of 
spectral index n = — 1.4 ± 0.1. 

Calculated at the lowest order, the values of cumulants of an unsmoothed field do not 
depend on the underlying power spectrum and are equal to the values calculated for a 
smoothed field with spectral index n = —3. This is not, however, the case for higher-order 
corrections to their values. In the case of scale-free power spectra (|76 ) when no smoothing 
is applied (W(kR) = 1) we obtain (see Appendix C) 

1 297 

£2 = ^ + h(n) « 0.3 (97) 

where h(n) is the part weakly dependent on the spectral index n which increases the 
rational number by roughly 10%. The remaining terms in the expression for the coefficient 
ax, equation (|93|), are skewness and kurtosis- related quantities and it was sufficient to 
calculate them at the lowest order. Combining equation fl9?D with the values of these terms 
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corresponding to the no smoothing case (i.e. the values for n = —3 in Table 2 and Table 5 
and the skewness values from equation (f79"l) ) we finally get 



Oi»l - 0.4 a 2 . 



(98) 



When smoothing is introduced, for n = — 



2 we obtain (see Appendix C) 




71 



0.369. 



(99) 



Combined with the other numbers calculated for n = —2 this result yields 



ai = 1-0.172 a 2 . 



(100) 



In the range —2 < n < — 1 we calculate £ 2 numerically and S 4 by interpolating the 
values given in Table 5. The calculation provides an independent check of the result for 
n = —2, equation (j99|), obtained analytically. Table 6 shows the two corrections to at 
separately: while the part containing S 4 remains roughly constant, £ 2 grows with n until 
it blows up at n = — 1. The last column of Table 6 lists the values of a\ in the cr-dependent 
way in order not to obscure the results by choosing arbitrary normalization needed for 
estimating a. Since the value of a 2 is of the order of unity on the scales of interest, it 
is clear from Table 5 that at the observationally preferred spectral index n « —1.4 the 
value of ai significantly departs from unity. It must be noted, however, that the nonlinear 
correction strongly depends on n and reaches zero between n = —1.6 and n = —1.7. 

To provide an example of the values of a\ we have normalized the power law spectra so 
that linear rms fluctuation in spheres of radius 8 h~ l Mpc is equal to unity. The resulting 
values of a 2 and a\ for spectral indices n = —1.4 ±0.1 at two different Gaussian smoothing 
scales R = 5 h- 1 Mpc and R = 12 hr 1 Mpc are listed in Table 7. 

Due to the reasons mentioned in the previous subsection we cannot explicitly examine 
the dependence of the coefficient a\ on Q. Still, it is constructed from moments involving 
second and third order solutions which have been proved to depend on Q very weakly. 
Consequently, the expected dependence of a\ on Q is weak. 

As an example of a scale-dependent power spectrum we again adopt the standard CDM 
model which, because of its behaviour at large wave-numbers {P{k) oc A: -3 ), does not 
introduce any challenges in the integration. The values of S 2 can be calculated numerically 
for a given smoothing scale. By combining with skewness values from Table 3 and the 
interpolated kurtosis-type values from Table 5 we end up with the coefficient a\ for the 
CDM spectrum. The values for different smoothing scales are given in the last column of 
Table 4. Although the S 2 values grow with scale (the remaining input to the correction to 
a\ remains roughly constant for this range of scales, see the last column of Table 5), the a 2 
values decrease much faster and, as we would expect for the perturbative results, at larger 
(i.e. more linear) scales the coefficient ai approaches its linear value, unity. 
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Tables 2, 4, 6 and 7 and Figures 1, 2 and 3 summarize the main results of this section. 
Table 2 provides the values of the coefficients a 2 and 03 for power law spectra in the whole 
range of the spectral index: —3 < n < 1. Those results are plotted in Figure 1 which also 
shows the ^-dependence of the coefficient a 2 . The values of the coefficient ai for power law 
spectra and the range of spectral index —2 < n < — 1 are given in Table 6. The correction 
to unity divided by a 2 is plotted in Figure 2. Table 7 lists the numerical values of a\ at 
n = — 1.4 ± 0.1 for two different smoothing scales when the <7g = 1 (top hat) normalization 
is adopted. The coefficients a±, a 2 and 03 for the standard CDM spectrum for a wide range 
of smoothing scales are provided in Table 4. Figure 3 shows their dependence on smoothing 
radius in the weakly nonlinear range of scales. 



4 Disentangling ft and linear bias 

The bottom line of our calculations is to propose a method, based on nonlinear corrections 
to the linear DVR, for measuring independently Q and bias from an IRAS-POTENT-like, 
density-velocity comparison. Let us assume that the galaxy and mass density contrast 
fields are linearly related, i.e. 5 g = b5, or 



We introduce a new variable 



By definition (Rf) we have 



5 = b- 1 5 g . (101) 



7^V-v. (102) 
-"0 



9 = f-\Q)5 v . (103) 



For the sake of simplicity let us consider the case in which a cosmic field is smoothed 
over a sufficiently large volume that the third order corrections to the weakly nonlinear 
DVR, equation (|65|), can be neglected. Using equations ( |101| ) and ( |103| ) we can rewrite 
DVR in the form relating two observables: the galaxy density contrast, 5 g , and the (minus) 
divergence of the velocity field, 5 V . We have 

(S 9 )\ Sv = jS v + a 2 j 2 (5 2 v -a 2 v ). (104) 

In the previous section we showed that the coefficient a 2 practically does not depend 
on Q. One can thus propose the following method for disentangling the effects of Q and 
linear bias. First, as the output of POTENT take simply 5 V (i.e. without any corrections 
for nonlinear ity ) . Next, plot the diagram 5 g -5 v . Finally, fit to the points a second order 
polynomial, 

5 g = Cl 5 v + c 2 {6 2 - a 2 v ) . (105) 
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Comparing equation ( |105|) with ( |104| ) we see that the fitted coefficients C\ and c 2 are related 
to / and b by 

ci = j, (106) 

and 

b . . 

C2 = a 2 j^. (107) 

So far, only the linear coefficient c\ has been measured. The results are usually expressed 
in terms of the variable (Dekel et al. 1993) 

P = c^. (108) 

The difficulty that from the linear density-velocity comparison one can estimate only 
the ratio f(Q)/b, or b/f(Q), is sometimes called the 'SVbias degeneracy problem'. In our 
opinion, however, there is nothing 'degenerated' in the fact that one cannot infer the values 
of two variables from only one equation involving them. True degeneration would happen 
if in the formula for the coefficient c 2 the parameters b and / entered only as a ratio, e.g. as 
(b/ f) 2 . Clearly, it is not the case. We can therefore solve equations ( |1U6| )- ( [ITJ7D separately 
for b and /. The result for / is 

/(fi) = ajj- = aa/TV . (109) 
c 2 

Unfortunately, the assumption of purely linear bias can be seriously questioned. In- 
deed, the value of the IRAS skewness is merely a half of the predicted one if b = 1 and this 
discrepancy is commonly attributed not to linear bias but to the fact that the IRAS survey 
systematically underestimates the density of galaxies in the cores of rich clusters (Bouchet 
et al. 1993). It simply means that the IRAS galaxies are nonlinearly (anti)biased tracers 
of mass distribution. In general, there is no a priori reason for the assumption of linear 
bias, justified for small (linear) fluctuations, to hold also in the case of fluctuations that 
are weakly nonlinear. Therefore, a correct method for estimating Q from the comparison of 
the weakly nonlinear galaxy density with velocity fields should take into account the effects 
of nonlinear bias. Such a method has been invented by Bernardeau (private communica- 
tion) and will be presented in the follow-up paper (Bernardeau, Chodorowski & Lokas, in 
preparation). 



5 Summary and concluding remarks 

In the present paper we derived a weakly nonlinear relation between cosmic density and 
velocity fields. In linear theory, the mass density contrast, 5, and the velocity divergence, 
V • v, are in a given point linearly related. If the fields are nonlinear the density-velocity 
relation (DVR) is neither linear nor local. In weakly nonlinear regime, however, the spread 
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around the mean trend is so small that the conditional mean can serve as a very useful 
local approximation of a true nonlinear DVR. 

We computed mean 5 given 9 := — /(fi) _1 iy ~ 1 V • v, that is (5)\ e , up to third order in 
perturbation theory. According to our calculations, it is given by a third order polynomial 
in 9, (5) \e = a\9 + a 2 (9 2 — a 2 ) + a 3 # 3 . This formula constitutes therefore an extension of the 
formula (5)\g = 9 + a 2 {9 2 — a 2 ), found by Bernardeau (1992), which involved second order 
perturbative solutions. Third order perturbative corrections not only introduce the cubic 
term but cause the a\ coefficient to depart from unity as well, in contrast with the linear 
theory prediction. We computed the numerical values of the coefficients a p for power-law 
spectra, as well as for standard CDM. The coefficients obey the hierarchy <^ a 2 <^ a±, 
which means that the perturbative series converges very fast. 

The key point of the method for disentangling Q and bias lies in the fact that the 
coefficients a p are practically ^-independent. We have shown that the dependence of the 
coefficient a 2 on Q is extremely weak. We have also given some arguments for the assump- 
tion that this is also the case for the coefficients a\ and 03. The detailed analysis of this 
problem will be given elsewhere. 

Recently Ganon et al. (1996) have performed a set of N-body simulations for a CDM 
family of models in order to test different local approximations to DVR in weakly nonlinear 
regime. They tested, among others, an approximation of 5 as a third order polynomial in 
9. By visual inspection of the plots they provide one can see that the approximation works 
excellently for \S\ less than unity. It does not surprise us, since it is just the regime of 
applicability of perturbation theory. The values of the fitted parameters are in qualitative 
agreement with our perturbative calculations for a standard CDM: a\ is slightly greater 
than unity, a 2 ~ 0.3 and a 3 is equal to a few hundredths. 

In order to derive accurately the values of a p from N-body one should, however, treat 
properly the final velocity field, determined in a simulation only at a set of discrete points 
(final positions). A two-step smoothing procedure, commonly used, leads to rather sub- 
stantial discrepancies between N-body simulations and analytical perturbative calculations 
of higher-order reduced moments (Lokas et al. 1995). Recently, Bernardeau & van de Wey- 
gaert (1996) proposed a new method for accurate velocity statistics estimation, based on 
the use of the Voronoi and Delaunay tessellations (adapted for a top-hat window function, 
however). The method proved to recover the tails of the velocity divergence distribution 
very accurately. Since the coefficients a p are given by skewness and kurtosis-like quantities 
(see section 2), probing the tails of the density and velocity distribution, the application of 
the method is necessary to recover the accurate values of the coefficients from simulations. 
This will be the subject of the follow-up paper (Bernardeau, Chodorowski & Lokas, in 
preparation). 
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Appendix A 

All of the perturbative calculations in Section 3 are much simpler if they are performed in 
Fourier space. For the first order of the density contrast field we have 



ik-x 



-ik-x 



(110) 



fill) 



<Ji(M) = D{t) J d 3 x <5i(x)e 

and the inverse Fourier transform is 

<?i(x,t) = D(t)(27T)- 3 J d 3 k 5i(k)e 

For the calculations of the coefficients a n the second and third order solutions for the 
density contrast and velocity divergence are needed and we give them here in the Fourier 
representation (e.g. Goroff et al. 1986). For the density field we have 

D 2 



5 2 (k,t) = ^ J d 3 p J d 3 q 5 D (p + q-k)5 1 (p)5 1 (q)P^(p,q) (112) 
* 3 (M) = J?-Jd 3 pJ d 3 q J d 3 r fo(p + q + r-k) < 5 1 (p) < 5 1 (q) ( y 1 (r)PS ) (p,q,r). (113) 



(2k) 

The symmetrized kernels for the density field are of the form 

p 2 ?(p> q) = + q> P' q) 



(114) 



4 s) (p,q,r) 



A 5 [#(p + q + r,p)J(q + r,q,r) + 
+ #(p + q + r,q + r)L(q + r,q,r)] + 
+B S F(p + q + r, p, q + r)L(q + r, q, r) + 

/ p q \ / p -»• r \ 



(115) 



+ 



q — > r 
\ r^p / 



+ 



q -> p 

V r^q / 



where A s = 1/108 and B 5 
Makino et al. (1992) i.e. 



1/189. In the expression above the notation follows that of 



#(p,q) 

^(p + q,p,q) 



p ■ q 

q 2 

1 |p + q| 2 p • q 

2 



J(p + q,p,q) = 4<E^_ 7 £z£p.,, i:i 



£(p + q,p,q) 



p 2 q 2 
p 2 q 2 



p 2 q 2 

p 2 + q 2 
7 — ^-p • q + 6. 

p A q A 



(116) 
(117) 
(118) 
(119) 



The solutions for the second and third order of the velocity divergence, 9 2 and 9 3 , are 
of the same form as the density solutions except for the kernels P$ and P$ that must 
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be replaced by the corresponding kernels P 2 e an d P30 ■ The second order kernel for the 
velocity divergence is 

^ S) (p ) q) = ^(p + q>p>q)- (120) 

The third order kernel for the velocity divergence is obtained from the density kernel ( |1 1 5| ) 
by replacing the constants As and B$ with Aq = 1/252 and B e = 1/63 respectively. 



Appendix B 



The expression (of 02^2) can be calculated in the same way as similar expressions (dfO^) 
or (SfS^)- This is in fact an easier part of the kurtosis calculations as it involves only the 
second order perturbative solutions. We have 



(%s 2 e 2 ) 



1 



i967r 3 r 3 (^ 



d 3 p / d 3 g / d 3 r P(p)P(q)P(r) x 



x W(|p + q|)W(|r - q\)W{p)W{r) x 
x J(p + q,p,q)L(r-q,r,-q) 



'12V 



which after performing the integrations over angular variables becomes 



a' 



r 3 i 



J dpjdq I dr (prr+^V+V^ 2 - 2 



x 



x 



26 , v ( q r\ 16 _ , ' 

— li(qr) — — | — -fa (or) H J 5 (or) 

21 2 W ; \r q 2 W ; 21 . 



Expanding the Bessel functions in powers 



Liz) 



E 

m=0 



m\T(y + m+ 1) V2 



i/+2m 



and using the fact that 



p x e~P dp=-T 



1 / 1 + sc 



(122) 



(123) 



io A \ A J 

we obtain the result as a series of Gamma functions which can be summed numerically up 

to arbitrary accuracy. The numerical results are given in the third column of Table 5. 



Appendix C 

Using the second order solution for the density field (p. 12 ) and corresponding one for the 
velocity divergence we obtain 



a 4 98(2vr)V 



J d 3 p J d 3 q P{p)P{q)W 2 {\p + q\R) M(p + q, p, q) (125) 
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where a and (3 stand for S or 6. For different combinations we have 



M 




for a 
for a 
for a 



(3 = 5 
5, (3 = 6 
(3 = 6 



(126) 



with J given by equation ( |118| ) and L by equation QH9| ). 
The second type of terms involves the third order solution 



(axa 3 ) 



J d 3 p J d 3 q P(p)P(q)W 2 (qR) 

{ A [ H{q,-p) J(p + q,p,q) 

+F(q,p + q) L(p + q,p,q)] 
- B F(q, -p, p + q) L(p + q, p, q)} . 



(127) 



If a = 5 the constants Ag and B$ must be used while if a = 6 they should be replaced 
respectively with A$ and Bq. The numerical values of the constants were given after 
equations (|115| ) and (|120|) respectively. 

After integration the expressions can be rewritten in a general way 



a' 



2ttV 4 



where 

P22(k) 

with 



k 3 



J™ dxP(kx) J +1 dfiP (k^l + x 2 - 2xfij — 



dk k 2 W 2 {kR)P ij {k) 



98(2vr) 



f(x,fx) 



2x/j,y 



(3x + 7/i- lOx/2 2 ) 2 
(3x + 7/i — 10x/x 2 )(7/x 
(7/i — x — 6x/j, 2 ) 2 



x 



6x/i 2 



for a 
for a 
for a 



(3 = 5 
5, (3 = 6 
(3 = 6 



and 



with 



^13 (AO 



k 3 p(k) 



dxP{kx)g(x) 



128) 



(129) 



(130) 



(131) 



9\P) 



i 

504 
1 

168 



l+x 



P - 158 + 100x 2 - 42x 4 + Jr(x 2 - l) 3 (7x 2 + 2) In ^ a 
i| - 82 + 4x 2 - 6x 4 + ^(x 2 - l) 3 (x 2 + 2) In p[±J 



for a 
for a 



= 5 
(3 = 6. 
(132) 

In the case of power law spectra, all the expressions ( |128j ) diverge individually in the 
limit of k — > and kx = q — » if n < — 1. Fortunately, as it is clear from the definition fl94}) , 
in calculating S 2 similar expressions are subtracted and all such diverging terms cancel out. 
In the opposite limit of k — > oo, q — > oo, divergencies occur for the terms involving second 
order if n > | and for the terms containing third order if n > — 1. As thoroughly discussed 
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by Lokas et al. (1996) the only way to get along with those divergencies at n > — 1 is to 
introduce a cutoff in the initial power spectrum at large wave-numbers. The results then 
depend on the cutoff wave-number k c . 

For integer values of the spectral index the integrals ( |128| ) can be performed by first 
finding the closed form expressions for similar to those proposed by Makino et al. 
(1992). In the limit of large cutoff wave-number k c — > oo a simple analytical result can 
be found by identifying the terms that dominate in this limit and integrating term by 
term. Then the result for n = —2 does not depend on the cutoff. 
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t H t {v) 

1 

1 v 

2 u 2 -l 

3 v 3 - 3v 

4 z/ 4 - 6z/ 2 + 3 

5 v 5 - 10z/ 3 + Ihv 

6 z/ 6 - 15z/ 4 + 45z/ 2 - 15 



Table 1: The Hermite polynomials 



spectral index n 


a 2 


a 3 


-3.0 
-2.5 


± w 0.190 
0.192 


3969 

-0.00935 


-2.0 


0.196 


-0.00548 


-1.5 


0.203 


-0.000127 


-1.0 


0.213 


0.00713 


-0.5 


0.227 


0.0165 





0.246 


0.0279 


0.5 


0.270 


0.0408 


1.0 


0.301 


0.0532 



Table 2: The coefficients a 2 and a 3 as functions of the spectral index n for scale-free power 
spectra and Gaussian smoothing 
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R 


U eff 


Ss8 




a 2 


a 2 (n eff ) 


5 


-0.946 


3.459 


2.177 


0.2136 


0.2141 


10 


-0.523 


3.305 


1.953 


0.2253 


0.2262 


15 


-0.255 


3.227 


1.821 


0.2343 


0.2355 


20 


-0.0654 


3.179 


1.729 


0.2416 


0.2430 


50 


0.462 


3.082 


1.482 


0.2666 


0.2680 


100 


0.722 


3.051 


1.359 


0.2819 


0.2830 



Table 3: The comparison of the values of the coefficient a 2 for the CDM spectrum calculated 
using the exact (fifth column) and the approximate (sixth column) method 



R 


a 2 


ai 


a 2 


&3 


5 


0.578 


1.196 


0.214 


0.00804 


10 


0.121 


1.119 


0.225 


0.0160 


15 


0.0419 


1.0755 


0.234 


0.0218 


20 


0.0185 


1.0519 


0.242 


0.0263 


50 


0.000975 


1.0125 


0.267 


0.0398 


100 


0.0000803 


1.00357 


0.282 


0.0465 



Table 4: The coefficients a±, a 2 and 03 for the CDM spectrum. In addition, the second 
column provides the values of the linear variance of the density (velocity divergence) field 
smoothed with a Gaussian filter for the CDM spectrum normalized as described in the text. 



spectral index n 




{5l5 2 6 2 )/a* 


(<WAr 6 




£ 4 


-3.0 
-2.5 


2.39 


If -4.01 
3.22 


1-3.61 
2.73 


f - 2.25 
1.53 


5536 
1323 

3.68 


-2.0 


1.87 


2.61 


2.11 


1.01 


3.31 


-1.5 


1.47 


2.13 


1.68 


0.631 


3.04 


-1.0 


1.16 


1.76 


1.38 


0.332 


2.84 


-0.5 


0.929 


1.47 


1.17 


0.0799 


2.71 





0.755 


1.24 


1.02 


-0.155 


2.63 


0.5 


0.638 


1.06 


0.919 


-0.398 


2.58 


1.0 


0.584 


0.916 


0.855 


-0.677 


2.53 



Table 5: The kurtosis-type quantities needed for the calculation of the coefficients a x and 
a 3 as functions of the spectral index n for scale-free power spectra and Gaussian smoothing 
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spectral index n 


s 2 


AS 31 W3 - £ 4 /2 


ai 


-2.0 


0.369 


-0.541 


1-0.172 a 2 


-1.9 


0.399 


-0.532 


1-0.134 g 2 


-1.8 


0.440 


-0.525 


1 - 0.0850 g 2 


-1.7 


0.496 


-0.518 


1 -0.0217 a 2 


-1.6 


0.576 


-0.511 


1 + 0.0643 o 2 


-1.5 


0.693 


-0.506 


1 + 0.187 g 2 


-1.4 


0.877 


-0.501 


1 + 0.376 g 2 


-1.3 


1.19 


-0.497 


1 + 0.698 g 2 


-1.2 


1.85 


-0.493 


1 + 1.36 g 2 


-1.1 


3.86 


-0.490 


1 + 3.37 a 2 



Table 6: The contributions to the weakly nonlinear correction to the coefficient ai in the 
most interesting range of spectral indices for power law spectra and Gaussian smoothing 



spectral index n R g 2 a\ 



-1.5 


5 


0.656 


1.12 


-1.4 


5 


0.638 


1.24 


-1.3 


5 


0.620 


1.43 


-1.5 


12 


0.176 


1.03 


-1.4 


12 


0.157 


1.06 


-1.3 


12 


0.140 


1.10 



Table 7: The values of the coefficient ai for the observationally preferred range of spectral 
indices of power law spectra and two different smoothing scales 
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Figure 1: The coefficients a 2 and 03 for scale-free power spectra and Gaussian smoothing 
as functions of the spectral index n. The solid lines correspond to the case of Q = 1. The 
coefficient a 2 is also shown for two other values of Q parameter (dashed lines). 
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Figure 2: The weakly nonlinear correction to the coefficient a\ divided by the linear vari- 
ance a 2 for scale-free power spectra and Gaussian smoothing in the observationally most 
interesting range of the spectral index n. 
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Figure 3: The coefficients a±, a 2 and a 3 for the standard CDM spectrum normalized to 
(top-hat) erg = 1 in the weakly nonlinear range of Gaussian smoothing scales. 
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